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Abstract This paper presents a study of the Poincare-Hough model of rotation of the 
synchronous natural satellites, in which these bodies are assumed to be composed of a 
rigid mantle and a triaxial cavity filled with inviscid fluid of constant uniform density 
and vorticity. In considering an Io-like body on a low eccentricity orbit, we describe 
the different possible behaviors of the system, depending on the size, polar flattening 
and shape of the core. 

We use for that the numerical tool. We propagate numerically the Hamilton equa- 
tions of the systems, before expressing the resulting variables under a quasi-periodic 
representation. This expression is obtained numerically by frequency analysis. This al- 
lows us to characterise the equilibria of the system, and to distinguish the causes of 
their time variations. 

We show that, even without orbital eccentricity, the system can have complex 
behaviors, in particular when the core is highly flattened. In such a case, the polar 
motion is forced by several degrees and longitudinal librations appear. This is due to 
splitting of the equilibrium position of the polar motion. We also get a shift of the 
obliquity when the polar flattening of the core is small. 

Keywords Natural satellites • Rotation • Periodic Orbits • Hamiltonian Systems • 
Numerical Methods 



1 Introduction 

Space missions like Galileo for the Jovian system or Cassini for the Saturnian one 
give us information on the internal structure of the natural satellites, through their 
gravity fields (Anderson 2001 pQ), observations of their surfaces (Porco et al. 2006 
39 ) or measurements of their rotational states (Tiscareno et al. 2009 |46| . Lorenz et 
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al. 2008 [28] ). It is known that the internal structure of a body influences its rotational 
dynamics, especially when this body is locked in a spin-orbit resonance, like the 1:1 
resonance for most of the natural satellites of the Solar system, and the 3:2 resonance 
for Mercury. 

There are at least two ways to approach the modelisation of the interactions be- 
tween the internal structure and the rotational dynamics. One way is to complexify 
the internal structure, taking account for instance of an atmosphere, a deformable 
crust, a subsurface ocean, an iron core. . . in a simplified dynamical model that allows 
to consider only one degree of freedom (see e.g. Rambaux et al. (2011) |41j for the 
longitudinal libration of satellites having an internal ocean, or Tokano et al. (2011) [UJ 
for the forcing of the polar motion of Titan due to its atmosphere). Another possibility 
is to consider a simple internal structure model (i.e. to assume the body to be rigid), 
in a full dynamical model considering several degrees of freedom (longitudinal motion, 
obliquity, and polar motion) like in (Henrard 2005 [19II20| ). 

An evolution of this approach is to consider a two-layer body composed of a rigid 
mantle and an ellipsoidal fluid core in which the flow is laminar and core-mantle interac- 
tions result in pressure coupling at the core-mantle boundary. This has been originally 
written by Hough (1895) [22] and Poincare (1910) 38 (that is the reason why this 
model is sometimes called the Poincare-Hough model), put in Hamiltonian form by 
Getino & Ferrandiz (see e.g. |13II15] ) under general assumptions, and by Touma & 
Wisdom (2001) [48], and recently used for Io (Henrard 2008 [2T]1. Mercury (Noyelles 
et al. 2010 [35]) and the Moon (Meyer & Wisdom 2011 [30]). Another model exists, 
taking account of the elasticity of the mantle (Getino & Ferrandiz 1995 Q3]). This case 
will not be considered here. 

In the case of the 1:1 spin-orbit resonance, the existing studies do not consider a 
wide range of internal structure parameter. This paper aims at contributing to fill the 
gap to understand the behavior of the system for any size and shape (provided it is 
triaxial) of the core. The plan of the study is the following: after a description of the 
model, we present a systematic numerical study of the system with different sizes and 
shapes of the core, the considered body being an Io-like body on a low eccentricity orbit 
with a uniform nodal regression and constant inclination. Then, "unusual" behaviors 
are highlighted with analytical explanations. 

2 The model 

In the study of Henrard |21] . the size of the core was not constrained, but its shape 
was assumed to be proportional to the whole Io. We here generalize this approach, in 
letting the shape parameters vary. 

2.1 Physical model 

Four references frames are considered (see Fig[T]&[2]). The first one, (ei,e2,e3) is 
assumed to be inertial for the rotational dynamics, it is in fact centered on the satellite 
and in translation with the inertial reference frame in which the motion of the satellite 
is defined. The second one, (n^, ri2, n§) is linked to the angular momentum of a pseudo- 
core that we define later, while the third one, i.e. (111,112,113), is linked to the total 
angular momentum of the satellite. Finally, the last one, written as (fi, f2, fs), is rigidly 
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Fig. 1 In the upper panel we have 3 reference frames: one linked to the ecliptic plane 
(ei,e2,e3), another linked to the angular momentum N (111,112,113), and the last one linked 
to the axes of inertia (fi, f2, fa) of the satellite. In the lower panel we have a similar configu- 
ration but instead of the angular momentum of the satellite, we have a reference frame linked 
to the angular momentum of a pseudo-core (defined later). We have the Euler angles (h, K, g) 
positioning the vector 112 on the plane perpendicular to the angular momentum of the satellite 
and the Euler angles (h c , K c , g c ) positioning the vector on the plane perpendicular to the 
angular momentum of the pseudo-core. The angles (I, J) and (l c , J c ) position the axis of least 
inertia. Note that J c is defined on the other side than J. 



linked to the principal axes of inertia of the satellite. In this last reference frame, the 
matrix of inertia of the satellite reads: 



(1) 




with < A < B < C, while that of the core is: 



(2) 

in the same reference frame. So, the orientations of the mantle and the cavity are the 
same, a misalignment of their principal axes would require to consider the mantle as 
elastic, this is beyond the scope of the paper. This would in fact require additional 
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parameters related to the elasticity of the mantle, see e.g. (Getino & Ferrandiz 1995 

El)- 

As for the whole satellite, we have < A c < B c < C c - In this way, the principal 
moments of inertia of the mantle are respectively A m — A — A c , B m — B — B c and 
Cm — C — C c - The principal elliptical radii of the cavity are written respectively a, b, 
c, yielding 

^(6 2 + c 2 ), 
^(a 2 + c 2 ), 
^(a 2 + 6 2 ), 

where p and M c are respectively the mass density and the mass of the fluid core, the 
quadrature being performed over the volume of the core. 



A, 



= JTJ (x\ + x\)pdx\ dx 2 dx 3 = 



B c = 
C c = 



(xi + x^)pdx\ dx2 dx% = ■ 
JT(a; 2 + x^£)pdx\ dx2 dx3 = ■ 



2.2 The kinetic energy of the system 

A Hamiltonian formulation of such a problem is usually composed of a kinetic energy 
and a disturbing potential, here the perturbation of the planet. Therefore, we consider 
every internal process, as the core-mantle interactions in our case, as part of the kinetic 
energy of the satellite. This section is widely inspired from (Henrard 2008 1211). 

The components (v\, 1)2,1)3) of the velocity field at the location (xi, X2, £3) inside 
the liquid core, in the frame of the principal axes of inertia of the mantle, are assumed 
to be (Poincare 1910 [55]): 



«1 = (w2 + -^2^x1 - (u 3 + x 2 , (3) 

v 2 = (u)3 + ~ u sjxi - (u>i + ^vij x 3 , (4) 
V3 = ^1 + j ) v l)x 2 - [tjj% + -U2J xi, (5) 

where (u)\, 0J2, ^3) are the components of the angular velocity of the mantle with respect 
to an inertial frame, and the vector of coordinates (y\, 1*2, V3) specifies the velocity field 
of the core with respect to the moving mantle. This vector is the velocity of a given 
fluid particle. Here we assume that this velocity field (i^i, ^2, ^3) depends only on the 
time t, and not on the spatial coordinates (asi, X2, X3). It implies that we have 

V v — ® V + ® V + ^ V — (6) 
dx\ dx2 8x3 

this equation is known as the continuity equation. 

The angular momentum of the core N^. is obtained by: 

Nc = JJJ( X x v)pdxi dx2 dx3 (7) 

core 

and the result is: 
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N' = 



+ 



Mc 
5 

Mc 
5 

M c 



C \ 2 / b \ 2 

-V\ + Ul I + I -V\ + UJl ] c 



C \ 2 / O \ 2 

-!/ 2 +^2j« +l-^2+^2)C 



6 \ 2 /(l \ ,2 

-ira + UJ 3 I a + I -I/ 3 + W 3 ) ft 



fl 
f 2 
f 3 - 



(8) 



We now set the following quantities: 



A = 



2M C 



6c = ,/(A c -B c + C c )(A c + B c -C c 



D 2 = HUac = ^/( - Ac + B c + C c ) {Ac + B C - C c ), 
-Ac + B c + C c ) (Ac ~B C + Cc), 



D 3 = HUab 



that have the dimension of moments of inertia and can be seen as parameters of the 
core as A c , B c and C c , and we can write: 



Nc = [Acuji + £>iz/i]fi + [B c u)2 + 2 ^]f2 + [C"c^ 3 + D 3 v 3 ]f 3 , 
while the angular momentum of the mantle is 

N m = A m wifi + B m i02h + C m u 3 h, 
and the total angular momentum of the satellite is 

N = [Awi + Divi] fi + [Blu 2 + D2V2] f 2 + [Ccj 3 + D 3 v 3 ] f 3 . 
The kinetic energy of the core is 



T c = i JJJ pv 2 dx\ dx2 dx 3 



(9) 



(10) 



(11) 



(12) 



T c = - (a c (Jx + uf ) + 5c(^l + + C c (cul + f 3 ) + 2D\w\v\ + 21)2^2^2 + 21)3^3^3 



while the kinetic energy of the mantle T m is 



T11 



1 Am^i + B m ^2 ~l~ CmOJ 3 

-N m • w . 



(13) 
(14) 



From T — T m + T c we finally deduce the kinetic energy of the satellite: 



T= -(Aui+Bi^2+CLji+Aciyi+Bciyi+Cc^+2D 1 i J j 1 u 1 +2D 2 i J j 2 u2+2D 3 uj 3 u 3 ). (15) 
We can easily check the expressions of the partial derivatives, for instance 



1 wc here correct a misprint present in Eq.13 of (Noyelles et al. 2010 35 ) 
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=Aui + D in =JVi (16) 



AT 

— = D 1 wi + J 4ci^=JVf, (17) 
wi 

where iVj are the components of the total angular momentum. N? are not the com- 
ponents of the angular momentum of the core but are close to it for a cavity close to 
spherical. We have, for instance for the first component: 

iVf - iV( c = (A c - Di)( Wl - I*) = ^(c - b) 2 K - u x ), (18) 

so the difference is of the second order in departure from the sphericity. From now on, 
we call angular momentum of the pseudo-core the vector N c = iVf fi + N^2 + iVf f3. 

With these notations, the Poincare-Hough's equations of motion, for the system 
mantle-core in the absence of external torque, are (see e.g. Eq.15 in [48] or \21\): 



^ = N x V N T, (19) 
^- = N c x V_ N cT, (20) 



with 



and 



v - r =S fi+ £ f2+ S f3 ' ^ 

v - NcT = ~ ~ ^k f3 ' (22) 

Here T is the kinetic energy expressed in terms of the components of the vectors N 
and N c , i.e. 



T = —(AcNt + A(Ni) ~ 2D x N x Nf) + — (B C 7V| + B(N%f ~ 2D 2 iV 2 7V 2 c ) 

2Q " ~ (23) 



+— {C c Ni + C(Ni) J - 2D 3 N 3 N§) 



with a = AA C - D\ , /3 = BB C - D\ and 7 = CC C - B\ 



2.3 The Hamiltonian 

2.3.1 The rotational kinetic energy 

We assume that the cavity and the satellite are almost spherical, this allows us to 
introduce the four small parameters eji 
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£1 
£2 
£3 
£4 



2C-A-B T MR 2 
2C = J2 



c 



2C 
2C C -A 



= 2C22 



Mir 

c 



Be 



2C C 
Be- A c 
2C C ' 



(24) 
(25) 
(26) 
(27) 



where M is the mass of our body and R its mean radius, and also the parameter 
S — C c /C, i.e. the ratio between the polar inertial momentum of the core and of the 
satellite. e\ represents the polar flattening of the satellite, while £2 is its equatorial 
ellipticity. £3 and £4 have the same meaning for the cavity. If we assume the core of 
the satellite to be spherical, we should take £3 = £4 = 0, while £4 = represents an 
axisymmetric cavity. Henrard [21 . considered that the ellipsoid of inertia of the core and 
the mantle were proportional, the mathematical formulation was £3 = e\ and £4 = £2. 

We now introduce the two sets of Andoyer's variables [2], (I, g,h, L, G, H) and 
(l c , gc, h c , L c , G c , H c ), related respectively to the whole satellite and to its core. The 
angles (h,K,g) are the Euler angles of the vector 112, node of the equatorial plane 
over the plane perpendicular to the angular momentum N, the angles (J, I) position 
the axis of least inertia fi with respect to 112. Correspondingly the angles (h c , K c , gc) 
are the Euler angles of the vector n^, node of the equatorial plane over the plane per- 
pendicular to the angular momentum of the pseudo-core N c , and (J c ,lc) position the 
axis of least inertia with respect to n§. The Figure [2] shows a schematic view of all the 
reference frames and relevant angles. The variables are (h,g,l) and (hc,g c ,lc) and the 
corresponding momenta (H — NcosK, G = N , L = TV cos J) and (H c — N c cos K c , 
G c = N c , L c = N c cos J c ). Expressed in Andoyer's variables the components of N and 
N c are: 



Ni = VG 2 -L 2 sinZ, 
N 2 = VG 2 - L 2 cos/, 
N 3 = L, 



JVf = VG 2 - L 2 sink, 
Ni = <s/G 2 - L 2 cosZc, 
NI, = L c . 



We can now straightforwardly derive the Hamiltonian Hi of the free rotation of the 
satellite, using Andoyer's variables and changing the sign of N c to take the minus sign of 
the Poincare-Hough equations into account (Eq |20|) . We also linearize the Hamiltonian 
with respect to the small parameters £^ (their orders of magnitude being about 10~ 5 , 
see Tab[T]), and get: 



s 




Fig. 2 The four reference frames gathered in the same view. The angles (h, K) position the 
plane orthogonal to the angular momentum N. The Euler angles (g, J, I) locate the axis of 
least inertia and the body frame (fi,f2,f3). The angles (J C ,Z C ) place the angular momentum 
of the pseudo-core with respect to the axis of least inertia f\ . 



1 . Gr 



Ho = 2C{1 _ 5) (G + -f + 2\/(G 2 - L 2 )(G 2 _ L 2 } cos( ; _ lc) + 2LLc 



c-l 



2C(l-8) 2 



G 2 -L 2 + G 2 C - L 2 + 2 J(G 2 - L 2 )(G 2 C - L%) cos(l - l c ) 



' 1 ' (G 2 - L 2 ) cos(2Z) + (G 2 - L 2 ) cos(22 c ) 



2C(1 -<5) 2 
+2^(G 2 -£ 2 )(G 2 -L 2 )cos(Z + ^ 

( 6(G 2 - L 2 ) + (G 2 C - L 2 C )(2 - i) 



2C(1 -8) 2 V ; v c cjy 8 

+28^/(G 2 -L 2 )(G 2 -L 2 )co S (l - l c )\ 



+ 2C(1~S) 2 ~ ^ COS(2 ° + (G ' ~ L ' )(2 ~ \ ] COS{2lc) 



+2Sy / (G 2 -L 2 )(G 2 -L 2 )cos(l + l c )\ . 



(28) 



We now introduce the following canonical change of variables, of multiplier n 
being the mean orbital motion of the satellite: 
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P = I + g + h, 
r — —h, 

Ci = - x /2P(l-cos J) sinZ, 
Pc = -Z c + 0c + h c , 

r c = -he, 

£2 = \/2P c (l + cos J c ) sin Z c , 



■R = P( l-cosJr), 

771 = y / 2P(l — cos J) cos Z, 

P c = P c (l- cos K c ), 

r/2 = \/2P c (l + cos J c ) cos/ c . 



The first three lines of this new set of variables and associated moments are related 
to the whole body, while the last three ones are related to the pseudo-core. P is the 
normalized norm of the angular momentum, it should be close to 1 at the spin-orbit 
resonance. Since the obliquity K is small, we have R oc K 2 , i.e. this is a small quantity 
related to the obliquity of the body. The quantities are related to the polar 

motion of the body, i.e. the angle J between the geometrical polar axis and the angular 
momentum, while I is the precession angle associated. We can note that £1 and r\\ are 
always defined, while Z is not defined when J = 0. The last three lines have basically the 
same meaning for the pseudo-core. We will see later that the degree of freedom (r c , R c ) 
is in fact not involved in the dynamics of this model, and that p c is not involved either, 
letting the norm of the angular momentum of the pseudo-core P c to be a constant. So, 
we can consider that the rotational dynamics of our body has 4, and not 6, degrees of 
freedom. 

In order to be consistent with the minus sign in the equations and before Z c , the 
amplitude of the wobble of the pseudo-core J c has to be replaced by tt — J c . In this 
way, we have L c — G c cos(7r — J c ) — — G c cos(J c ). In this new set of variables, we have 



N^nC^pi-jp-^) 2 ^, 
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and the Hamiltonian of the free rotational motion becomes, after division by nC: 



«(p-a±i)(a±i-^)) 

+ ^ - ^) (a -2±i)(^ -«.&>) 
.^(p-2+3i)( R -a+i)( M+{l6 )) 

+ 2S<I (P - S±5) (ft _ („, 2 + £l£2 ) j . (30) 



Finally, in order to get an easy-to-read formula, we can develop this Hamiltonian 
up to the second order in (£i, £2, »7l, to get: 
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+ 2(1 "^ )2 ( p(Zi + m) + Pc(d + m) + 2Vpp c ( mm - 
+ 2(1 "1 2 J)2 ( p (?i " + Pc fe 2 - %) - 2^PP c ( mm + ) (31) 



2 



2(l-5) 2 



+ 2 



<">(£ 2 + ^ 2 ) + (2 - J) + ^) + 25VPPc{vim - 



This is in fact a third-order development since the powers in £2, ?7l, ^2) are even. In 
the forthcoming computations, this last approximation has not been used, the equations 
we have propagated deriving from the Hamiltonian (|30p . 

2.3.2 The gravitational potential 

To compute the gravitational potential due to the parent planet on its satellite, we must 
first obtain the coordinates x, y, and z of the unit vector pointing to the planet in the 
reference frame linked to the principal axes of inertia (fi,f2,f3), from its coordinates 
in the inertial frame x^, yi and 2j. Five rotations are to be performed: 

(yy^j = R 3 (-l)Ri(-J)R3(-g)Ri(-K)R 3 (-h) (j^j (32) 

with Xj, yi, Zi depending on the mean longitude A , the longitude of the ascending 
node fl a , the longitude of the perihelion zu , the inclination i, and the eccentricity e. 
The rotation matrices are defined by 

/ cos </> - sin cf> \ / 1 \ 

Rai'P) ~ I s' m 4> cos(f> , Rl(4>) ~ I cos0 — s'm(f> . (33) 
\ 1 / \ sin 4> cos <t> ) 

The gravitational potential then reads: 

V X {X , 1, g, h, J, K) = -\C Q -^- (6i (x 2 + y 2 ) + e 2 (x 2 - y 2 )) (34) 

where Q is the gravitational constant, M p the mass of the perturber, i.e. Jupiter for 
Io, [x, y, z) the unit vector pointing at the perturber in the frame (fi, f%, fa), such that 
x 2 + y 2 + z 2 = 1, while d is the distance planet-satellite. 

Let us note that unlike Henrard [21], we consider that the perturbation is applied to 
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the whole satellite and not only to its mantle, this issue is addressed in Noyelles et al. 
(2010 [35]). 

From the variables x, y and z, it is easy to introduce the set of variables defined in (Eq. 
1291) . We also modify the moment A associated with A (that appears in the expressions 
of x and y) in such way that all our variables are now canonical with multiplier 1/nC 
and our gravitational potential becomes (after division by nC) 

1-i2{\o 1 p,P,r,R,ii,ri 1 ) = ~~pF( ei ( x2 + y 2 "> + _ f 2 ))' ( 35 ) 

Finally, we use the formulae ()30|) and (I35|) to get the Hamiltonian of the system: 

H = Hi (P, ft , m , ft, m) + H 2 (A , p, P, r, R, ft , 77!). (36) 

The four degrees of freedom of this Hamiltonian are the spin (p, P), the obliquity (r, 
R), the wobble of the whole body (ft, 771) and the wobble of the core (ft, 772)- 

2.3.3 Evaluating P c 

Since the variable p c , spin angle of the pseudo-core, does not appear explicitly in 
the Hamiltonian of the system, its associated momentum P c , norm of the angular 
momentum of the pseudo-core is not ruled by the Hamilton equations. So, it can be 
either a constant, or a time varying input as is the orbital motion of the system. We 
here choose to set P c — S = C c /C, the mean value of P being very close to 1 as our 
pseudo-Io is in 1:1 spin-orbit resonance. So, we assume a kind of equipartition of the 
norm of the angular momentum between the core and the mantle. 

An exact equipartition would be Pc(t) = SP(t), meaning that the fluid would follow 
every fluctuation of the orbital velocity of our pseudo-Io. It would mean that the fluid 
follows the longitudinal librations of the mantle, as if it were rigid. In such a case, the 
amplitude of the longitudinal librations would not be affected by the presence of an 
at least partially liquid core. Observations of such librations for Mercury (Margot et 
al. 2007 [12]) and the Moon (Koziel 1967 [M|, Williams et al. 1973 [SO]) support the 
assumption that the longitudinal librations are the response of the solid mantle (and 
not of the full body) to variations of the orbital velocity of the body. That is the reason 
why we consider a constant value for P c , that results from a kind of rough averaging 
of P. 

While this model describes the rigid dynamics of a body having a fluid core, we 
must not forget that real bodies on which this model could be applied have a viscous 
fluid core. We here discuss the relevance of our assumptions on P c for these bodies. 
From a physical point of view, the reason for the decoupling between the fluid and the 
mantle is a low viscosity of the fluid. At the core- mantle boundary (CMB), the no-slip 
condition imposes that the velocity field follows the mantle. So, there is a thin turbulent 
layer close to this boundary, known as the Ekman layer, in which the velocity field 
evolves continuously from the no-slip condition at the boundary to the one satisfying 
P c = 8. The typical thickness of the Ekman layer is d = \Jv / Q (Greenspan 1968 |17|). 
v being the kinematic viscosity and O — n the spin frequency of the fluid. Usually a 
kinematic viscosity v = 10~ 6 m 2 /s is considered at the core-mantle boundary because 
it is consistent with a Fe/Fe-S composition (see e.g. Kerswell 1998 [23]), what yields 
d — 0.16 m. A viscosity of 36m 2 /s is necessary for the thickness of the Ekman layer 
to reach 1 km. In fact, the viscosity is expected to increase with the depth under the 
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CMB, since molten, and even rigid iron, should be concentrated at the inner core (see 
e.g. Rutter et al. |44|). We anecdotally recall the extremum of viscosity of the pitch 
derived from the pitch drop experiment set up in 1927 at the University of Queensland, 
Australia (Edgeworth et al. 1984 QJ]), i.e. v = (2.09 x 10 5 ±4.6 x 10 4 )m 2 /s. 



2.4 Link with the Navier-Stokes equation 

As said in the introduction, there are at least two ways to approach the interactions 
between the internal structure and the rotational dynamics. One is to complexify the 
internal structure in considering only one degree of freedom, and the other one is to 
consider several dynamical degrees of freedom (4 in this study) with a quite simple 
internal structure. We must keep in mind that these 2 very different approaches aim 
at studying the same bodies. A complete study of the core dynamics would require to 
consider the Navier-Stokes equation, we here make a link with this equation to help in 
the interpretation of our model from a physical point of view. 

The dynamics of a particle of fluid is often assumed to be ruled by the well-known 
Navier-Stokes equation, we give here its expression as given in (Greenspan 1968 [17]): 

-|q+(q-V)q+2fixq=-ivp-i/Vx (Vxq)-rx^, (37) 

with 

— q: particle velocity measured in a rotating system 

— f2: angular velocity of the rotating system, its coordinates being (u)i, wq, W3) 

— p: density of the fluid 

— r: position of the particle 

— p = P + pU — ^(f2 x r) ■ (17 x r): the reduced pressure, where P is the pressure of 
the fluid, and U an exterior potential, 

— v. kinematic viscosity of the fluid. 

In our case we have 

/ (a/c)v2X3 - (0/6)1/3x2 \ 
q = (b/a)u 3 xi - (b/c)vix 3 . (38) 
\ (c/b)v x x2 - {c/a)v 2 x 1 J 

In an over-simplified case where we neglect the viscosity v, the convective acceler- 
ation (q • V)q and the reduced pressure p, the formula (|37[) reads: 

^q + 2fixq = 0, (39) 



+ 2(oj 2 ^3 - = 0, 

^ + 2(^3^ - Wl i>3) = 0, (40) 
+ 2((JiV2 - W2"l) = 0. 
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For comparison, the formula ()19[) reads: 




(Bu)2 + Div-i)uij, 



(Aji + D\Vi)ljJ2 



(Cuj 3 + D 3 u 3 )uji 



(Blo 2 + D2V2)U\. 



(Aoji + Diui)oj 3 , 



(Cu 3 + D 3 v 3 )u 2 , 



(41) 



The systems of equations (|40p and (|4ip present some similarities, the main difference 
being that the moments of inertia are involved in Eg 1411 They should be in fact con- 
sidered as global equations (i.e. considering the whole volume of fluid), while the Eg 1401 
is a local form, considering an individual fluid particle. 

The reader can find another formulation of these equations in (Rambaux et al. 2007 



3 A numerical study 

3.1 The algorithm 

As shown in Henrard [21], the proper frequency associated with the core, i.e. the free 
core nutation, is close to the spin period of the considered body. For a synchronous 
satellite, this period is also the orbital period, so we have a proximity between a 
proper frequency of the problem and a forcing period. As a consequence, a perturbative 
approach will meet difficulties to converge because of small divisors. Such a problem 
has already been encountered in (Noyelles et al. 2010 [35])- That is the reason why 
we prefer a full numerical study, consisting of a numerical integration of the equations 
derived from the Hamiltonian (|36p . and a frequency analysis of the solutions of the 
problem. The frequency analysis algorithm we use is widely inspired from NAFF (see 
Laskar 1993 |26| for the method, and Laskar 2005 [27] for the convergence proofs), with 
a refinement suggested by Champenois (1998 0) consisting in iterating the process to 
enhance the accuracy of the determination. 

The basic idea of the frequency analysis is to consider that a complex variable of 
the problem x(t) is quasi-periodic, i.e. can be expressed as a, a priori infinite, sum of 
a converging trigonometric series like 



n=0 

where A n are constant complex amplitudes, and v n constant frequencies, with 



TO- 



oo 




(42) 



N 




(43) 



the bullet meaning that the coefficients have been numerically determined. A detailed 
description of the algorithm is given in appendix. In the case of a real variable, the 
Eg 1431 becomes 
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x 



N 

if) ~ A * n cos (y' t + 0* ) > 

n=0 



(44) 



or 



x 



N 

(t) ~ A ' n sin + ' 

71=0 



(45) 



where the amplitudes are now real, and the <j>n are rea l phases expressed with the 
counterclockwise convention, previously included in the complex amplitudes. 

The rotation of a synchronous satellite is reputed to have reached an equilibrium 
state, known as Cassini State 1 (see e.g. Cassini 1693 [3], Peale 1969 [37], and Bouquillon 
et al. 2003 [3] for an extension to the polar motion), after dissipation of its rotational 
energy. There should remain free oscillations with negligible amplitude around the 
equilibrium, in the following we assume them as null, since they cannot be detected 
except for the Moon (Rambaux & Williams 2011 [12]). It can be shown that, for rigid 
dynamics, between 2 and 4 Cassini States exist. In the context of natural satellites of 
the giant planets where the nodal precession rate is small with respect to the orbital 
frequency, the 4 Cassini States exist, and they induce an obliquity close to fcS, k 
being an integer (see Ward & Hamilton 2004 g5] or Noyelles 2010 [33], Appendix B). 
The Cassini State 1, corresponding to k — 0, i.e. a small obliquity, is a priori the most 
probable one, because it is stable and the primordial obliquity of the satellite is thought 
to be small. 

In order to numerically simulate the rotational dynamics of the satellite, we need 
initial conditions that are actually very close to the equilibrium, that is perturbed by 
the orbital dynamics of the satellite. For that, we use the algorithm NAFFO (Noyelles 
et al. 2011 [35]), consisting in: 

1. A first numerical integration of the equations of the system, with initial conditions 
conveniently chosen, 

2. Frequency analysis of the solution and identification of the contributions depending 
on the free modes, 

3. Evaluation of the free modes at the time origin of the numerical simulation, and 
removal from the initial conditions, 

then the process is iterated until convergence. In a Hamiltonian framework as is the 
case here, Noyelles et al. [35] have shown that the convergence is quadratic in the 
amplitude of the free modes, provided that the quasi-periodic decomposition is exact, 
i.e. that the signal is indeed quasi-periodic, and that the numerical error has negligible 
impact. The proof is based on the d'Alembert characteristic (see e.g. Henrard 1974 
18 ), that gives a relation between the amplitudes A n and the arguments u n t in Eq |43l 
This algorithm has already been successfully applied in problem of rotational dynamics 
(Dufey et al. 2009 [11], Noyelles 2009 [33], Robutel et al. 2011 [43]), in dynamics of 
exoplanetary systems (Couetdic et al. 2010 [8]), and in the analysis of ground-track 
resonances around Vesta (Delsate 2011 [5]). 

3.2 The numerical tests 



The numerical algorithm we have just described has been used in different cases, de- 
pendent on the free parameters £3, 64 (« polar flattening and equatorial ellipticity of 
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Table 1 Physical and dynamical parameters ruling our pseudo-Io. We used the same as 21 . 
The orbital frequency n and the regression rate of the ascending orbital node £2 are taken 
from LI. 2 ephemerides | 25l . The phases (initial conditions of the orbital angles) are arbitrarily 
chosen. 



Parameter 


Value 


GM P (planet) 


1.261648547674763616 x 10 23 fcm 3 .s~* 


GM (satellite) 


5955.5 fcm 3 .s -2 


R p 


71492 km 


J2p 


1.4736 x 1CT 2 


7o 

'J2 


1 S2S x 1 o~ 3 


C*22 


5.537 x 1(T 4 


C/(MR 2 ) 


0.376856 


a 


422029.958 km 


e 


4.15 x 10" 3 


I 


2.16 arcmin 


n 


1297.2044725279755 rad/y 


w 


0.97311853791375 rad/y 


tl 


-0.8455888497945 rad/y 


A o (0) 





U7 o (0) 


2 rad 


n o (0) 


0.1 rad 


t MR 2 

ei = J2—Q— 


4.85066 x 10~ 3 


or , MR 2 
e 2 = 2022—77— 


2.93852 x 10~ 3 



the core), and S = C c /C, representing the size of the core through its inertial polar 
momentum. In all our simulations we considered a kind of pseudo-Io, i.e. a satellite 
with physical and dynamical properties close to the ones of the Galilean satellite of 
Jupiter Io, except that its orbit has constant eccentricity and inclination. The numerical 
integrations are performed with the Adams-Bashforth-Moulton 10th order predictor- 
corrector integrator, with a tolerance of 10 _1 , and a step size of 5x 10~ 5 y w 1.8 x 10 — 3 
d. 

We considered as reference values for the internal structure parameters 8 = 0.5, 
£3 = ei and €4 — €2, and we tested different pseudo-Ios with different values of these 
parameters. 



4 "Classical" behavior 

We expect to have, at the Cassini State 1: 

— a = p — X + tt close to because of the 1:1 spin-orbit resonance, 

— P close to 1 (the norm of the angular momentum being close to nC), 

— p — £l a — h — S7 D + r (third Cassini Law), 

— R close to (the obliquity being small), 

— J and J c close to (small polar motions of the satellite and its core), 

the "classical" behavior being small oscillations around this equilibrium. We use it to 
define our first initial conditions, before refining them with NAFFO. 
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Table 2 Proper frequencies of the small oscillations around the equilibrium for £3 = ei, 
£2 = £4 and S = 0.5. n is the orbital frequency given in Tab[l] 





Frequency 


Period 


ui/n 




(rad/y) 


(d) 






243.4050908 


9.42845 


0.187638 


LU V 


4.1898509 


547.73630 


3.2299 x 10- 3 


Ld w 


19.5319416 


117.49643 


0.015057 


Ul z 


1334.4264821 


1.71979 


1.028694 



Table 3 The variable P-l. The series are in cosine. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 






1 


1.5156914 x 10" 4 


1296.2313540 


65.408° 


1.77047 


A — TO 


+ 7T 


2 


6.6760683 x 10~ 7 


2592.4627080 


-49.183° 


0.88523 


2A — 2to 


+ 7T 


3 


3.8653845 x 10" 9 


3888.6940620 


-163.775° 


0.59016 


3A D — 3to 


+ 7T 


4 


1.2702319 x 10~ 9 





-180° 


00 




est 


5 


2.2957822 x 10" 11 


5184.9254160 


81.634° 


0.44262 


4A - 4to d 


+ 7T 


6 


1.5952860 x lO -11 


2596.1001228 


-11.459° 


0.88399 


2A Q - 


2fi 



4.1 In-depth study of a reference case 

We here present an in-depth study of a "reference case", considering £3 = ex, £4 = £2> 
and 8 — 0.5. This study consists of a numerical estimation of the frequencies of the 
proper librations (Tab[2]), and of a numerical decomposition of the canonical variables 
(TabOtoil). 

We recall that the orbital frequency n is 1297.20447137 rad/y (Lainey et al. 2006 
25 ). A comparison with Henrard 21j lacks of significance since the physical model was 
different (the gravitational torque of Jupiter acting only on the mantle, while it acts on 
the whole satellite here), but we can see that, like Henrard, we find a proper frequency 
of the core u z close to the spin frequency of Io, that is also its orbital frequency since 
our satellite is locked in the 1:1 spin-orbit resonance. 

The Tab[3]to[5]give a quasi-periodic decomposition of the canonical variables with 
identification of the forced oscillations, i.e. the mean longitude of our pseudo-Io A , 
the motion of its pericenter zo , and the motion of its orbital ascending node £l a . The 
phases are indicated at the time origin and allow to determine the presence of n or 7r/2 
in the identification. 

Since our rotational model takes 4 degrees of freedom into account, we can split 
the canonical variables and moments into 3 groups related to these degrees of freedom. 

The first group (a, P) (Tab[3] & [4]) can be linked to the longitudinal motion. We 
can see that the mean position is the theoretical equilibrium (a — 0, P — 1) related 
to the 1:1 spin-orbit resonance, and there are small oscillations around this equilib- 
rium, related to the mean anomaly A — tjJo and its harmonics. We can see from the 
Tab[4]that the deviation from the theoretical equilibrium does not exceed 2 arcmin for 
an eccentricity of 4.15 x 10 -3 . This amplitude is proportional to the eccentricity (at 
least for small eccentricities, see e.g. Comstock & Bills 2003 [7]), that induces periodic 
variations of the planet-satellite distance. 

The second group (p,R) (Tab[5] & [6]) locates the angular momentum of the whole 
body with respect to the orbital plane. Once more, the angle can be averaged to 
with a instantaneous departure that does not exceed 2 arcmin, this equilibrium is a 
consequence of the third Cassini law. It is also known that the mean obliquity, that can 
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Table 4 The resonant argument a. The series are in cosine. 





Amplitude 


Frequency 


Phase 


T 




Identification 






(rad/y) 


(t=0) 


(d) 








1 


62.574 arcsec 


1296.2313540 


-24.592° 


1.77047 


A 


o — Too 


+ 7T/2 


2 


0.138 arcsec 


2592.4627080 


-139.183° 


0.88523 


2A 


— 2to 


+ 7T/2 


3 


0.053 arcsec 


3888.6940620 


106.225° 


0.59016 


3A D 


— 3to 


+ 7T/2 


4 


2.03 X 10~ 5 arcsec 


2596.1001228 


-101.459° 


0.88399 


2A D 


-2fi 


- tt/2 


5 


2.37 X 10~ 6 arcsec 


5184.9254164 


-8.366° 


0.44262 


4A D 


— 4ro 


+ 7T/2 


6 


1.21 X 10~ 6 arcsec 


1299.8687682 


166.868° 


1.76551 


A + Too 


-2fi D 


- tt/2 



Table 5 The variable R. The series are in cosine. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 




1 


2.5966515 x 10" Y 





0° 


CO 


est 


2 


1.6920424 x lO" 10 


2596.1001228 


-11.459° 


0.88399 


2A - 2«o 


3 


3.1006440 x lO" 11 


1296.2313540 


65.408° 


1.77047 


A — Too + 7T 


4 


1.2914058 x 10~ 12 


3892.3314767 


-36.051° 


0.58960 


3A D - Too - 2^ + t/2 


5 


5.5914871 x 10" 13 


3.6374149 


-142.276° 


630.924 


2to - 2«o 



Table 6 The variable p. The series are in cosine. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 




1 


67.204 arcsec 


2596.1001228 


78.541° 


0.88399 


2A - 2a o + t/2 


2 


1.542 arcsec 


1296.2313540 


155.408° 


1.77047 


Ao — Too — 7r/2 


3 


0.517 arcsec 


3892.3314767 


53.949° 


0.58960 


3A — Too — 2J7 + 7r 


4 


0.222 arcsec 


3.6374146 


-52.276° 


630.924 


2to d - 2«o - 3tt/2 



Table 7 The variable rji + The series are in complex exponential. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 








1 


5.22646 x 10~° 


-1298.0500614 


-174.270° 


1.76799 


-Ao 


+ «o 


— 7T 


2 


5.54231 x 10" 7 


1298.0500614 


174.270° 


1.76799 


A 


-Slo 


+ 7T 


3 


2.81640 x 10" 7 


-1.8187074 


71.138° 


1261.849 




— Too 


+ 7T 


4 


2.20597 x 10~ 7 


1.8187074 


-71.138° 


1261.849 




-fio 


— 7T 


5 


5.96421 x 10" 9 


-2594.2814154 


-59.679° 


0.88461 


Too + £lo 


- 2A 


— 7T 


6 


3.66084 x 10" 9 


2594.2814154 


59.679° 


0.88461 


2A D — Too 


-Slo 


+ 7T 



be derived from the mean value of R, is due to the interior structure and the regression 
rate of the orbital node (see e.g. Ward & Hamilton 2004 [49]). We can also see that 
the oscillations are dominated by the mode 2A D — 2J2 OI emphasizing an influence of the 
orbital node on this degree of freedom. 

The third group involves the last two degrees of freedom 771) (Tab[7J) and (£2, V2) 
(Tab[SJ), that are strongly coupled as shown by Henrard |21j). They represent respec- 
tively the polar motion of the whole body and the orientation of the velocity field of 
the fluid. They are ruled by two kinds of small oscillations: fast ones due to harmonics 
of the proper mode A c — £l a , and slow ones due to the argument of the pericenter 



19 



Table 8 The variable r\2 + i(,2- The series are in complex exponential. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 




1 


7.36054 x 10"° 


1298.0500614 


-5.730° 


1.76799 


— S7 


2 


3.33184 x 10" 7 


-1298.0500614 


5.730° 


1.76799 


-Ao + tto 


3 


2.03442 x 10~ 7 


1.8187074 


108.862° 


1261.849 


- S7o 


4 


1.55904 x 10" 7 


-1.8187074 


-108.862° 


1261.849 


£l - ro 



Table 9 Influence of the size of the core <5, with £3 = £1 and £4 = £2. 



<5 


T u (d) 


T v (d) 


T w (d) 


T* (d) 




R* 


0.1 


12.650 


453.259 


208.551 


1.745 


2.304 x 


10" 


/ 


0.2 


11.926 


480.369 


185.790 


1.741 


2.365 x 


10- 


7 


0.3 


11.156 


504.669 


163.028 


1.736 


2.450 x 


10- 


7 


0.4 


10.328 


526.965 


140.264 


1.729 


2.525 x 


10- 


7 


0.5 


9.428 


547.734 


117.496 


1.720 


2.597 x 


10- 


7 


0.6 


8.433 


567.278 


94.723 


1.706 


2.680 x 


10- 


7 


0.7 


7.303 


585.780 


71.939 


1.685 


2.760 x 


10- 


7 


0.8 


5.963 


603.294 


49.130 


1.645 


2.842 x 


10- 


7 


0.9 


4.216 


619.394 


26.230 


1.540 


2.922 X 


10- 


7 



Table 10 Influence of the polar flattening of the core £3, with 8 = 0.5 and £4 = £2- 



£3/ei 


T u (d) 


T v (d) 


T w (d) 


T z 


(d) 








R* 


0.2 


9.428 


6414.819 


117.118 


1 


728 


1 


040 


X 


10" b 


0.3 


9.428 


2491.673 


117.112 


1 


727 


4 


612 


X 


10- 7 


0.4 


9.428 


1572.784 


117.121 


1 


726 


3 


592 


X 


10- 7 


0.5 


9.428 


1163.452 


117.147 


1 


725 


3 


177 


X 


10- 7 


1 


9.428 


547.734 


117.496 


1 


720 


2 


597 


X 


10- 7 


3 


9.428 


254.876 


122.879 


1 


695 


2 


326 


X 


10- 7 


5 


9.428 


210.742 


136.657 


1 


668 


2 


277 


X 


10- 7 


6 


9.428 


200.875 


148.926 


1 


655 


2 


268 


X 


10- 7 


7 


9.428 


194.454 


168.492 


1 


641 


2 


260 


X 


10- 7 


8 


9.428 


189.284 


201.639 


1 


628 


2 


256 


X 


10- 7 


9 


9.428 


185.621 


278.943 


1 


616 


2 


251 


X 


10- 7 



4.2 Influence of the parameters 



To characterise the influence of the internal structure parameters (i.e. £3, £4 and 8), 
we quantify their effects on our outputs. We choose here to consider in particular the 
proper frequencies lj u to tv z , and the mean value of R (Tabl9lto [TT1) . 

We can see that all these outputs depend on the size of the core 5 (Tab(5|. In 
particular, the period of the free longitudinal librations T u follows the classical law 
(see e.g. Goldreich & Peale 1966 [16]): 



(46) 



yielding T u oc y/1 — 8. We note that this period depends on the size of the core, while 
Henrard did not find any dependency in applying the gravitational torque just on the 
mantle. To check the influence of the shape of the core, we now present the outputs 
with varying e 3 (Tab[l0]) and e 4 (Tab[TT). 

We recall that for Mercury, i.e. in the case of the 3:2 spin-orbit resonance, the 
flattening of the core £3 alters the frequencies u) v and u z , but not the others ones. 
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Table 11 Influence of the equatorial cllipticity of the core £4, with <5 = 0.5 and £3 = ei. 





T u (d) 


T v (d) 


T w (d) 


Tz (d) 


R* 





9.428 


545.949 


117.771 


1.7199 


2.5996 x 10" 


7 


0.1 


9.428 


546.128 


117.718 


1.7199 


2.5998 x 10" 


7 


0.5 


9.428 


546.841 


117.563 


1.7199 


2.5954 x 10" 


7 


1 


9.428 


547.734 


117.496 


1.7198 


2.5967 x 10- 


7 


3 


9.428 


551.316 


118.652 


1.7195 


2.6070 x 10- 


7 


5 


9.428 


554.914 


122.283 


1.7193 


2.6069 x 10- 


7 


10 


9.428 


564.010 


149.248 


1.7186 


2.6248 x 10- 


7 



Here, the variations of the period of the free longitudinal librations T u have only 
negligible variations, while the 3 other proper frequencies are affected. As for Mercury, 
the periods T v and T z increase with £3 getting closer to 0, T z getting closer to the 
spin period 1.769 d, and T v tending to infinity. We also have an increase of the free 
wobble period when £3 increases. We can note that it seems to be possible to fine-tune 
the parameters (£3 ~ 7.7ei) to have a resonance between the free wobble and the free 
oscillations of the obliquity (T v — T w ), but this is only anecdotal. This very peculiar 
case would require strict fine-tuning between the flattening of the body and of the core 
to occur, so we can consider it as very unlikely. Finally, the equilibrium position of the 
angular momentum, i.e. R* , is shifted from the origin (here the normal to the orbit) 
when the core tends to be spherical (small £3). 

In the case of the 3:2 spin-orbit resonance, no significant influence of the equatorial 
ellipticity of the core had been detected. We here (Tab llip see a small influence on T v , 
T w , T z and R* , but that does not seem to be significant. Once more, the longitudinal 
librations seem not to be affected. 

As a reminder, Henrard ,19, found free periods of respectively T u = 13.25, T v — 
159.39 and T w = 229.85 days in considering a rigid Io. The rigid value of 13.25 days 
can be obtained in setting C m = C in Eg 1461 

5 Analysis of a bifurcation 

In the previous section, we do not present the behavior of the system for some physically 
possible values of the core shape parameters 63 and £4. The reason is that for some 
range of these parameters, the system presents more complex behaviors, that we here 
introduce. In particular, we assume since the beginning a "classical" Cassini State 
1 in which null amplitudes of the polar motions of the body J and of the core J c 
define a stable equilibrium. In fact, this has not been checked yet, and our numerical 
investigations have revealed that this equilibrium is unstable for instance for 8 = 0.5, 
£3 = 10ei and £4 = 0. 

5.1 Numerical characterisation of the equilibria 

A simulation of the behavior of the system for S = 0.5, £3 = lOfi and £4 = gives a 
butterfly shape for the outputs related to the polar motion of the body (jji, £1) and of 
the velocity field of the fluid (772 > £2) for the solution passing close to the equilibrium 
defined by J = J c = (Fig(3]). It suggests that this equilibrium is in fact unstable, and 
that two new stable equilibria appear. 
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Fig. 3 Trajectory passing close to the equilibrium £1 = £2 = »?i = V2 for 8 = 0.5, £3 = 10ei 
and £4 = 0. The left panel shows the polar motion of the whole body, and the right one is related 
to the pseudo-core. We can see that the trajectory does not librate around this equilibrium, 
but presents a butterfly-shape, that suggests the presence of 2 new stable equilibria. 



Table 12 The variable P-l for £3 = 10ei and £4 = 0. The series are in cosine. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 








1 


2.81032 X W-'-s 








OO 






est 


2 


1.56939 x 1(T 4 


1296.2313540 


65.408° 


1.77047 


A 


— TO D 


+ 7T 


3 


6.77223 x 10~ 7 


2592.4627080 


-49.183° 


0.88523 


2A - 


- 2to 


+ 7T 


4 


5.20072 x 1CT 8 


1298.0500614 


174.270° 


1.76799 


A 


-flo 


+ 7T 


5 


3.80792 x 10~ 8 


1.8187074 


288.862° 


1261.849 


TO 


-Slo 


+ 7T 


6 


3.92165 x 10~ 9 


3888.6940620 


-163.775° 


0.59016 


3A - 


- 3zu Q 


+ 7T 


7 


3.50511 x lO" 10 


2594.2814154 


59.679° 


0.88461 


2A Q — to 


-Slo 


+ 7T 



Table 13 The resonant argument a for €3 = 10ei and £4 = 0. The series are in cosine. 





Amplitude 


Frequency 
(rad/y) 


Phase 
(t=0) 


T 

(d) 




Identification 


1 


63.973 arcscc 


1296.2313540 


-24.592° 


1.77047 




Ao 


- Too 


+ 7T/2 


2 


0.139 arcsec 


2592.4627080 


-139.183° 


0.88523 




2A - 


2to 


+ 7T/2 


3 


0.113 arcsec 


1298.0500614 


84.270° 


1.76799 




Ao 


-Slo 


+ 7T/2 


4 


5.4 X 10 -4 arcsec 


3888.6940620 


106.225° 


0.59016 




3A D - 


3zu 


+ 7T/2 


5 


3.0 x 10 -5 arcsec 


2594.2814154 


149.679° 


0.88461 


2A 


— Too — 


«o + 


3tt/2 


6 


2.2 X 10 -5 arcsec 


1294.4126466 


46.546° 


1.77295 


Ao - 


2to + 


Ro + 


3tt/2 


7 


1.6 X 10 -5 arcsec 


2596.1001228 


-101.459° 


0.88399 




2A - 


2fi 


- tt/2 



These equilibria have been reached thanks to NAFFO. The quasi-periodic decom- 
positions of the solution corresponding to the equilibrium (£1 = £2 = 0, 771 « 0.25, 772 ~ 
—0.17) are given in Tab ll2l to ll7l The other equilibrium is symmetrical to this one, i.e. 
corresponds to (£1 = £2 — 0, rji ~ —0.25, 772 ~ 0.17). 

We can see from these tables that the difference is not only in (£1, r]\, (3, 772). The 
difference for the degree of freedom related to the longitudinal behavior (a, P) is strik- 
ing. First, we can see a significant departure (2.81 x 10 -3 ) from the expected mean P, 
i.e. 1 (Tab ll2p . We also note significant longitudinal librations related to the combina- 
tion of proper modes Ao — £l (Tab ll3p . that did not appear in the "classical" behavior 
(TabBl). 

The difference is even more important for the degree of freedom related to the 
location of the angular momentum, i.e. (p,R) (Tab ll5l fc I14[l . In this case, we can see 
large oscillations associated with the argument of the pericenter vj — fio- K i s known 
that a motion due to the position of the pericenter has the eccentricity as physical 
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Table 14 The variable R for €3 = 10ei and €4 = 0. The series are in cosine. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 




1 


3.3172485 x 10~ Y 





0° 


00 


est 


2 


3.0899861 x 1CT 7 


1.8187074 


108.862° 


1261.849 


TOo — fio 


3 


1.4728954 x lCT 10 


2596.1001228 


-11.459° 


0.88399 


2A D - 2« 


4 


7.9302319 x 10~ n 


2592.4627080 


-49.183° 


0.87986 


2A — 2to + 7T 


5 


2.0015617 x 10- 11 


1296.2313540 


65.408° 


1.77047 


A — TO + 7T 


6 


1.5218917 x 10- 11 


2594.2814154 


-120.321° 


0.88461 


2A — TO — fl a 


7 


8.1094453 x 10~ 12 


1294.4126466 


46.546° 


1.77295 


Xo - 2to d + fi D + 3tt/2 


8 


4.4761150 x 10" 12 


1298.0500614 


174.270° 


1.76799 


^0 " fit, + T 



Table 15 The variable p for €3 = 10ei and €4 = 0. The series are in cosine. 





Amplitude 


Frequency 
(rad/y) 


Phase 
(t=0) 


T 

(d) 






identification 


1 


39.163° 


1.8187074 


18.862° 


1261.849 




w 




-7T/2 


2 


13.384° 


3.6374148 


127.724° 


630.924 




2TO 


-2fi 


- ?r/2 


3 


6.099° 


5.4561222 


-123.414° 


420.616 




3to 


-m a 


- 7I-/2 


4 


3.127° 


7.2748300 


-14.552° 


315.462 




4ra D 




- tt/2 


5 


1.710° 


9.0935369 


94.310° 


252.370 




5to 




- ?r/2 


6 


58.431 arcmin 


10.9122443 


-156.828° 


210.308 






-6fi 


-tt/2 


7 


34.233 arcmin 


12.7309517 


-47.966° 


180.264 




7ra D 




-tt/2 


8 


20.474 arcmin 


14.5496590 


60.896° 


157.731 




8to 




- tt/2 


9 


12.440 arcmin 


16.3683661 


169.758° 


140.205 




9ra 


-m a 


- tt/2 


10 


7.653 arcmin 


18.1870733 


-81.380° 


126.185 




10ro o - 




- ?r/2 


11 


4.755 arcmin 


20.0057804 


27.482° 


114.714 




llro - 


-llfio 


- ?r/2 


12 


2.979 arcmin 


21.8244872 


136.344° 


105.154 




12ro - 


-12f2o 


- tt/2 


13 


1.880 arcmin 


23.6431936 


-114.794° 


97.065 




V&VJo - 


-130, 


- tt/2 


14 


1.715 arcmin 


2596.1001226 


78.541° 


0.88399 




2X0 


-2fi 


+ TT/2 


15 


1.193 arcmin 


25.4618992 


-5.932° 


90.132 




14ro - 


-140, 


- tt/2 


16 


1.168 arcmin 


2597.9188300 


-172.597° 


0.88337 


2A D 


+ to - 


-m ~ 


- 3tt/2 


17 


52.836 arcsec 


2594.2814149 


-30.321° 


0.88461 


2A 


— to 




- Ztt/2 


18 


47.694 arcsec 


2599.7375373 


-116.265° 


0.88276 


2A D - 


r 2to d - 


-4fi D - 


- 3tt/2 


19 


45.673 arcsec 


27.2806037 


102.930° 


84.123 




15to - 


-150, 


- tt/2 


20 


32.468 arcsec 


2601.5562447 


45.127° 


0.88214 


2A G - 


- 3TOo - 


-5fi - 


- Ztt/2 


21 


29.268 arcsec 


29.0993068 


-148.208° 


78.866 




16to - 


-16fi D 


- tt/2 


22 


22.103 arcsec 


2603.3749520 


153.989° 


0.88152 


2A D - 


r 4ro D - 


-6fi D - 


- 3tt/2 


23 


18.829 arcsec 


30.9180086 


-39.346° 


74.226 




17to - 


-17fi 


- tt/2 


24 


15.046 arcsec 


2605.1936596 


-97.149° 


0.88091 


2A D - 


- 5ro - 


-7f2 - 


- 3tt/2 



Table 16 The variable r/i + i£i for £3 = 10ei. The series are in complex exponential. 





Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 




1 


0.2501568 





0° 


00 


est 


2 


3.36571 x 10~ o 


-1296.2313540 


-155.408° 


1.77047 


to — X — 3tt/2 


3 


7.44538 x 10~ 6 


1296.2313540 


155.408° 


1.77047 


X — TO + 3-7T/2 


4 


5.19058 x 10- 6 


-1298.0500614 


-174.270° 


1.76799 


-Ao + fi - TT 


5 


1.72624 x 10~ 6 


1.8187074 


-71.138° 


1261.849 


TOo - fi D + TT 


6 


1.69988 x 10~ D 


-1.8187074 


71.138° 


1261.849 


fio - ro o - TT 


7 


5.90500 x 10~ 7 


1298.0500614 


174.270° 


1.76799 


Ao - fio + T 


8 


4.69722 x 10~ 8 


2592.4627080 


-49.183° 


0.88523 


2A — 2to + TT 


9 


3.95648 x 10~ 8 


-2592.4627080 


43.183° 


0.88523 


2to — 2A — TT 


10 


3.31650 x 10~ 9 


2594.2814161 


59.679° 


0.88461 


2X — TOo — fl + TT 
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rable 


17 The variable r\2 


+ i£,2 for £3 = 10ei 


The series 


are in complex 


exponential. 






Amplitude 


Frequency 


Phase 


T 


Identification 






(rad/y) 


(t=0) 


(d) 






1 


0.1690165 





180° 


CO 




est 


2 


3.25350 x 10" 5 


1296.2313540 - 


-114.592° 


1.77047 


A - 


TO 


3 


7.01986 x 10~ 6 


1298.0500614 


-5.730° 


1.76799 


A - 


fio 


4 


1.16880 x 10" 6 


-1.8187074 - 


-108.862° 


1261.848 


a D - 


Too 


5 


1.15291 x 10" 6 


1.8187074 


108.862° 


1261.848 


TO D - 


fto 


6 


1.09791 x 10" 6 


-1296.2313540 


-65.408° 


1.77047 


TO — A 


— 7T 


7 


2.37221 x 10" 7 


-1298.0500614 


5.730° 


1.76799 


flo" 


- A 


8 


8.29840 x 10~ 9 


2592.4627079 


130.817° 


0.88523 


2A D - 


2ro 


9 


4.26308 x 10" 9 


2594.2814150 


59.679° 


0.88461 2A 


o — Too - S7 


+ 7T 


10 


1.54598 x 10" 9 


-2592.4627077 - 


-130.817° 


0.88523 


2ro Q — 


2A 



cause, while our eccentricity is only 4.15 x 10 , the peak-to-peak oscillations of p 
reaching 80°. So, we can expect higher oscillations for bigger eccentricities. 

In this case, the shift of P led us to change iteratively the value of the constant 
P c so that it remains equal to 5 < P >. We have seen that a change of P c yields a 
significant difference on the locations of the stable equilibria, that is the reason why 
the mean values of r\\ + i£i and 772 + 1^2 we give in Tab ll6ll7l are significantly different 
from the ones that can be guessed from FigO 

5.2 Analytical study 

In order to understand the appearance of 2 new stable equilibria, we propose a simpli- 
fied analytical study of the problem. This study consists in starting from the Hamilto- 
nian H (Eq|36} , in expressing the oscillating angle (respectively a = p — X + tt because 
of the 1:1 spin-orbit resonance, and p — S7 — h because of the third Cassini Law), 
in averaging over the circulating ones, to deduce a secular Hamiltonian yielding the 
equilibria. All these calculations have been performed thanks to Maple software. 

The starting point is the Hamiltonian H (Eg 1361) in which the coordinates of the 
perturber (i.e. a pseudo- Jupiter if we consider a pseudo-Io) x and y are replaced thanks 
to Eg 1521 with 

Xi = - (cos fi cos(A - £l Q ) - cos I Q sin fi Q sin(A 
Vi = - (sinfi cos(A -fi D ) + cos I cosfi Q sin(A Q 
Zi = -sin7 sin(A - fi D ). 

We here neglect the influence of the eccentricity. 
Then the following canonical transformation is performed 

o = p — A + tt, P, 

p = Zlo + r , R, 
£1, m, 

Since this transformation, involving A D and £l a , is time-dependent, we must add —nP+ 
£IR to the Hamiltonian. a and p are oscillating arguments that can be averaged to 0, 
while A and £7 are circulating. 



■Slo)): 
■Slo)): 



(47) 
(48) 
(49) 



(50) 
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A first-order averaging of the Hamiltonian is performed, then the Hamilton equa- 
tions are derived, i.e. 



da an 

W — ~5P' 


dp an 
~m ~aiji 


dp an 
dt ~ 7m> 


dR au 

dt ~~ Tfr-i 


dti _ an 
~w ~~ cRjT' 


dm _ an 


d& _ an 

dt Wpi 


dm an 
dt - ~a~& 



(51) 



the equilibria corresponding to null time derivatives of the variables and associated 
moments, i.e. the right-hand side of these equations vanish. The numerical exploration 
drove us to neglect the influence of the inclination and the obliquity (1 = 0, R — 0), 
and to consider £i and £2 as null at the equilibrium. These approximations allowed 
us to simplify the system, and we finally find with a good agreement the equilibrium 
values of P, rji and 7/2 in solving numerically the following equations: 

1 da P — P c m , . . 77? .„„. 

n~dt = ~T~S~ + 2(1^? (£1 " 62 " *» + 5ei) + W=Sj (52) 

+ V1V2 / ej _ £2 _ Se3 + Se4 



2(1-5) ^PP C - P4/4 - Pctf/4 + 



1-5 



;t-(i^<«'--*' + '«» + ^ + 5i i^ 1 (-. + « + *.-*d- 5 ^5 

Cl — ^2 — Se 3 + Se4 

1-5 



4(1-5) ^JPPc - Pril/ 4 - P c rj{/ 4 + 7?^|/16 
Y^y/PPc - Pr)ll4 - P c rH/ A + V f V 2 2 /16 (l + 



and 



n dt (1-«5) 2 V \s -y-'v sjy'i-s 



+ w^w r 1 + £2 + 1 2 * ) e3 + U - 2 J 64 J - (54) 



0?i/ 4 - p ) A + ei - £2 -fe 3 + ^ 4 \ 

4(1 - S)yJ PPe - P7?|/4 - Peril! ± + VlVlf 16 ^ 1 " 6 ' 

+ Y=- 5 \fPPc ~ Prft/A ~ Peril/ 4 + tJtJ/1& (l + 



£1 — £2 ~ 3^3 + $£4 



1-8 

For £3 = 10ei, £4 = and 8 = 0.5, the real roots of this system are 

P = 1.046772470, r/i = 1.446908787, r? 2 = -1.023119016 
P = 1.002812138, m = 0.2502391659, r? 2 = -0.1690724173 
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Table 18 Location of a new stable equilibrium, determined analytically (a) thanks to Eg l52l 
to 1541 and numerically (n), for <5 = 0.5. The last column, A, gives the amplitude of the 1.77-d 
longitudinal librations, obtained in our numerical code. Here, only the equilibrium correspond- 
ing to i)i > and r?2 < has been considered. In all these cases, another stable equilibrium 
exists in changing the signs of tji and 772 ■ 



£3/ei 


£a/c2 


P - 1 (n) 




P - 1 (a) 




Vl (n) 


m ( a ) 


m (n) 


V2 (a) 


9.45 





1.5831 x 10" 


4 


4.0680 x 10" 


-4 


0.0608 


0.0972 


-0.0411 


-0.0657 


10 





2.8103 x 10" 


3 


2.8121 x 10" 


-3 


0.2502 


0.2502 


-0.1690 


-0.1691 


10 


0.3 


1.9482 x 10" 


3 


1.9501 x 10" 


-3 


0.2099 


0.2100 


-0.1418 


-0.1419 


- p = 


1, Vl = 


m = 

















- P = 1.002812138, 771 = -0.2502391659, 772 = 0.1690724173 

- P = 0.3489241565, 771 = 0.8353731582, 772 = -0.5606980247 

while they are, for £3 = 9ei, £4=0 and 8 = 0.5: 

- P = 1.041484268, 771 = -1.443249317, 772 = -1.020531379 

- P = 0.3471614224, 7/1 = -0.8332603709, 772 = -0.5892040580 

- P = 1, 771 = 772 = 

- P = 0.9978852209, 771 = -2.010693353, t? 2 = -1.421011855. 

So, we can see for P « 1 and \rji\, 1 772 1 < 0.5, an appearance of 2 additional equilibria. 
In order to test the validity of this analytical study, we propose fTab ll8[) a short com- 
parison between its results and the numerical results, in 3 cases where the 2 equilibria 
appear. We can see a significant discrepancy for the first case, where £3 = 9.45ei and 
£4 — 0. In this case, the equilibria are close to the origin 771 — 772 = 0, while a good 
agreement is reached for the other two cases, where the equilibrium values of 771 and 
772 are bigger. The observed discrepancy can be due to the neglect of the obliquity, the 
inclination and the eccentricity. 

We now propose to study the existence of these 2 additional equilibria. Since their 
existence is linked to the stability of the equilibrium corresponding to r]i, { j = 0, P = 1 
and P c = 8, we in fact study this stability. In setting fi = £2 — 0, P = 1 and P c — 8 
in the averaged Hamiltonian, we get the quantity S: 



5(771,772) = a - 1 + 



2(1-8) 



+ ei 




3 rfi + vlS - (vf + 77!) /4 
2 + 2(1 - 5) 2 




+ £2 



3 Vl + V26 - {nt + 4) /4 



(55) 



2 2(l-5) 2 



+ (^3 




-84 + 4(1 -28) + 84 /A 



) 



2(l-5) 2 



-8 



with 



mm . L 7725 + 772 t? 2 ??! 

= T^V — + ^6" 



(56) 
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We do not call S "Hamiltonian" since two variables, i.e. £1 and £2, are sets to constants, 
while their associated momenta r/i and 772 vary. The study is now equivalent to the 
investigation of the extrema of the surface defined by the Eq |55l In fact we study the 
point defined by 771 = 772 = 0, we know thanks to previous calculations that it gives null 
first-order derivatives of S. The topological nature of this point can be investigated in 
studying the second order partial derivatives of S. We consider the Hessian matrix 



M 




(57) 



(6- 



ei — ^2 + <5(e4 — £3 + 1 — 5) \/S(l — 5 + ei — e 2 + 5(e 4 — £3)) 
V5(l-8 + e 1 -e 2 + 5(e 4 -e 3 )) l-S + S(e 1 -e 2 -2e 3 + 2e 4 ) 



A minimum (corresponding to a stable equilibrium) is reached when the two eigenvalues 
of the Hessian, A1.2, are positive. We have: 



A 2 = /3- 



(58) 
(59) 



with 



1 - 5 2 + [e x - e 2 )(l + S) + (ea - e 4 )(l - 35) 



(60) 



and 



zi = - 5 2 ) - 2(ei - e 2 )(l -15 + 76 2 - 6 3 ) + 2(e 3 - e 4 )(l - 35 - 8 2 + 35 3 ) 

+ (ei - e 2 ) 2 (l + 5) 2 + 2(ei6 4 + £2^3 - ^3 - - 25 + 55 2 ) 

+ (e 3 -€ 4 ) 2 (l- 25 + 5 2 +45 3 ). (61) 



Numerical evaluations show that Ai is always positive, and that A2 is usually positive, 
except for the interior parameters given in Tab ll8l In these peculiar cases, we have 
A1A2 < 0, so the considered point (771 = 772 = 0) is a saddle point. 

This study shows that the equilibrium corresponding to J = J c — is unstable for 
A2 = 0. This condition is independent of the mean motion and is applicable to any 
body in 1:1 spin-orbit resonance, in which the interior model of a rigid mantle, a fluid 
core and a small solid inner core composed of dense material is realistic. We have also 
neglected the effect of the orbital inclination and of the regression of the ascending 
node. This approximation is relevant, since most of the natural satellites of the giant 
planets have inclinations of the order of a few arcmin, and the nodal regression of Io 
is one of the most rapid in the Solar System. Since here this approximation gives good 
results, it should be available for most of the Solar System bodies in a comparable 
dynamical situation. 
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Table 19 The first case is the reference one, and the other ones correspond to the cases where 
two additional stable equilibria appear. 



£3/^1 


e4/e2 


< K m > 


< Jm > 


< Jc > 


1 


1 


2.299 am 


0.155 as 


21.470 as 


9.45 





9.065 am 


3.632° 


3.335° 


10 





37.542 am 


14.975° 


13.726° 


10 


0.3 


31.220 am 


12.555° 


11.516° 



5.3 Effect on the observable variables 

We now consider the influence of this peculiar behavior on the observable parameters, 
i.e. data that could be observed if our pseudo-Io were real and if it were observed with 
enough accuracy. In particular, they have to refer to the mantle since its rotation is 
actually the rotation of the surface. These observable data can be deduced from the 
canonical variables, that give a complete mathematical description of the system. 

A complete derivation of the observable outputs can be found in (Noyelles et al. 
|35|). we here choose to represent the following quantities: 

— the mean obliquity of the mantle < K m >, 

— the mean amplitude of the polar motion of the mantle < J m >, 

— the mean amplitude of the polar motion of the core < J c >■ 

All these results are obtained thanks to frequency analysis, and they are gathered 
in Tab lliJI We can see that the stable equilibria that appear induce a forcing of the 
polar motion of the surface (or mantle) of our pseudo Io (Fig(4}, that can reach 15°. 
In (Noyelles 2008 [32]) we found a forcing of the polar motion of a rigid Titan, due to 
a resonance between the free wobble and the forced precession of Titan's perihelion. 
We considered it as a possible explanation for the super-synchronous rotation of Titan, 
before it was observed (Stiles et al. 2008 & 2010 [45])- This is different here, since no 
resonance appears. 




<5 = 0.5, e 3 /ei = t4/e 2 = 1 6 = 0.5, e 3 /ei = 10, e 4 = 

Fig. 4 Location of the North Pole of the mantle of the body (located by fa) with respect to 
its angular momentum N m in the classical case (left) and with a highly flattened core (right). 
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Table 20 The variable Kexp(ir) for <5 = 0.5, £3=0 and 62 = £4. The series are in complex 
exponential and the amplitudes in arcseconds. We can note a nearly constant component v, 
that has a negligible in the usual case. 





Amplitude 


Frequency 


Phase 


T 


Identification 




(arcsec) 


(rad/y) 


(t=0) 


(d) 




1 


84.516 


1.093 x 10" 4 


78.801° 


2.1 x 10' 


V 


2 


66.484 


0.8455888 


-5.730° 


2714.006 


-9,o 


3 


0.022 


-2595.2545339 


5.730° 


0.884 


Zlo - 2A 


4 


0.006 


1296.2314632 


-35.742° 


1.770 


A — mo + v 


5 


0.006 


-1296.2312447 


-166.558° 


1.770 


ro — A + v 


6 


0.005 


-1295.3857651 


-71.138° 


1.772 


ro — A — £7 D — T 


7 


0.005 


1297.0769428 


59.679° 


1.769 


A — VJ — £l a + TT 



6 Orientation of the angular momentum 

Among the Third Cassini Law (see e.g. Cassini (1693) 4 or Colombo (1966) [B]), the 
equilibrium orientation of the total angular momentum of the body is assumed to be 
in the Cassini State 1. As a consequence, the angular momentum, the normal to the 
orbital plane and the normal to the Laplace Plane are coplanar, the Laplace Plane 
being a reference plane based on the precessional motion of the orbital ascending node, 
that minimizes the variations of the inclination of the considered body. There are in 
fact several ways to define this plane, as for instance in (Yseboodt et al. 2006 [51]) or 
in (D'Hoedt et al. 10 ). A difficulty is: how to consider a constant reference plane if 
the precession rate of the ascending node is not constant? Should we average over a 
"long enough" time interval, or over a time-interval suitable to the observations of a 
space mission? 

The reader can find in (Noyelles 2009 [33] ) a discussion on the choice of an "ap- 
propriate" reference plane depending on the variations of the orbital inclination, that 
allows the argument p — £l D — h to librate. It is shown that, for the rotation of a rigid 
body in 1:1 spin-orbit resonance, if the satellite orbits close to its parent planet, the 
precessional motion is ruled by the oblateness of the planet (its J2) and so its preces- 
sion rate is close to be constant. In such a case, choosing the equatorial plane of the 
planet as a reference plane to describe the behavior of the angular momentum of the 
body can be a convenient choice. However, when the satellite orbits far from its parent 
planet as it is the case for Titan or Callisto, the reference plane for the nodal precession 
is shifted because of the Solar gravitational perturbation. In such a case, considering 
the planet's equatorial plane as the reference plane could either result in a oscillating 
rotation node h as it is the case for Titan (Noyelles et al. 2008 )31|). either result in an 
erratic apparent behavior due to an improper choice of the reference plane, as is the 
case for Callisto (Noyelles 2009 [55]). 

In our case of a pseudoTo with a constant regression of the node, no "strange" 
behavior is expected. In particular, the Tab[6j supports the assumption of a quasi- 
periodic behavior of the difference of the nodes p. However, we have found a different 
behavior for a small flattening of the core £3 (Fig[S]and Tab l20p resulting in a significant 
shift of the mean equilibrium orientation of the total angular momentum. This shift 
seems to be not constant but a long-period oscillation, the period being w 57, 000 years. 
We call v this oscillation. 

In (Noyelles et al. 2010 [35]), we had found a particular behavior for small £3, that 
we attributed to the exact resonance between the Free Core Nutation frequency uj z 
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-150 -100 -50 50 100 150 -60 -40 -20 20 40 60 80 100 

Kcos(r) (arcsecond) Kcos(r) farcsecond) 

5 = 0.5, £3 = ei, £4 = 62 8 = 0.5, £3 = 0, £4 = £2 

Fig. 5 Behavior of the orientation of the angular momentum of our pseudo Io K expir, with 
2 different internal structure models, in the inertial reference frame. The right panel shows a 
shift of this motion that is not on averaged at the origin. 




■150 -100 -50 50 100 150 -150 -100 -50 50 100 150 



KnjCosfr,,,) (arcsecond) KnjCosfr,,,) (arcsecond) 

S = 0.5, £3 = £1, £4 = £2 & = 0.5, £3 = 0, £4 = £2 

Fig. 6 Behavior of the orientation of the angular momentum of the mantle (i.e. the surface) of 
our pseudo Io -K"m expir m , with 2 different internal structure models, in the inertial reference 
frame. Contrary to the total angular momentum (Fig[5j, it does not exhibit particular behavior. 



and the spin frequency. We also noticed an asymptotic behavior of the free frequency 
ui v that tended to (and the period T v to infinity) when £3 tended to 0. This last 
behavior is here observed as well as can be seen in Tab llOl This is confirmed by some 
tests at £3 = ei/10 suggesting T v — 9933.75 days. However, even if the free period T z 
gets closer to the spin period of 1.76799 day, it does not seem to reach it. So we cannot 
speak of resonant behavior, it seems more likely to be a kind of singularity at £3 = 0. 

The Fig[S] shows the orientation of the angular momentum of the mantle/surface, 
that does not exhibit this shift. So, if such a situation would occur (i.e. very small polar 
flattening of the core), the equatorial/ring plane of the planet could be an acceptable 
reference plane to describe the orientation of this axis. In fact, a physical signature of 
this dynamics remains in the core, we indeed get a mean J c of « 3 arcmin for £3 = 
while we have < J c >~ 21 arcsec for £3 = ei. 



7 Conclusion 



In this study we have presented the behavior of a pseudo-Io orbit on a low eccentric orbit 
around its parent planet, with a uniform nodal regression and a constant inclination, 
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in considering it as a two-layer body composed of a rigid mantle and a fluid triaxial 
core. This model can be applied to study the rotation of most differentiated natural 
satellites. 

We have described the "usual" case, consisting of small oscillations around the 
expected equilibrium, i.e. synchronous rotation with a small obliquity and no polar 
motion, but we also have, especially for a highly flattened core, another behavior re- 
sulting in a polar motion forced by several degrees. Another peculiar behavior is when 
the polar flattening of the core is very small. In this last case we have a forcing of the 
obliquity of the full body, but not of its mantle, so there should be no observational 
evidence of this phenomenon. From a mathematical point of view, this could be due 
to a kind of singularity in the parameter 63. 

This study aimed at exploring the behavior of a model, its application to real 
bodies would require to consider complete ephemerides. This would add additional 
forcing frequencies complicating the dynamics of the system. New behavior cannot a 
priori be excluded. 

A possibility to improve the model would be to consider nonlinear phenomena in 
the fluid, but this is another story. . . 
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A The NAFF algorithm 

The frequency analysis algorithm that we use is based on Laskar's original idea, named NAFF 
as Numerical Analysis of the Fundamental Frequencies (see for instance Laskar 1993 26 for the 
method, and Laskar 2005 27] for the convergence proofs). It aims at identifying the coefficients 
ajz and ajj, of a complex signal f(t) obtained numerically over a finite time span [— T;T] and 
verifying 

n 

/(*) ~ ak ex P(*<^*)> ( 62 ) 
fe=l 

where oj^ are real frequencies and complex coefficients. If the signal f(t) is real, its frequency 
spectrum is symmetric and the complex amplitudes associated with the frequencies wj. and 
— LU^ arc complex conjugates. The frequencies and amplitudes associated are found with an 
iterative scheme. To determine the first frequency ui\ , one searches for the maximum of the 
amplitude of 

(j>(w) =< f(t), exp(«urt) >, (63) 
where the scalar product < f(t), g(t) > is defined by 

< /(*),*(*) >= ^ f f(t)g(t)* X (t)dt, (64) 

g(t)* being the complex conjugate of g(t). \(t) is a weight function alike a Hann or a Hamming 
window, i.e. a positive function verifying 

h f_ T x(t)dt = L (65) 

Using such a window can help the determination in reducing the amplitude of secondary 
minima in the transform 1164 l l. Its use is optional. 
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Once the first periodic term exp(iajit) is found, its complex amplitude a\ is obtained 
by orthogonal projection, and the process is started again on the remainder fi(t) = f(t) — 
ai exp(iu>it). The algorithm stops when two detected frequencies are too close to each other, 
what alters their determinations, or when the number of detected terms reaches a limit set by 
the user. This algorithm is very efficient, except when two frequencies are too close to each 
other. In that case, the algorithm is not confident in its accuracy and stops. When the difference 
between two frequencies is larger than twice the frequency associated with the length of the 
total time interval, the determination of each fundamental frequency is not perturbed by the 
other ones. Although the iterative method suggested by Champenois Jj allows to reduce this 
distance, some troubles may remain. In our particular case, these problems are likely to arise 
because of the proximity between the free frequency of the core ui z and the frequency of the 
spin. 
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